Week 3: Basic Geospatial Operations

PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2024

Author

Affiliation

Jeff Jacobs

jj1088@georgetown.edu

Open slides in new tab →

How to Do Things with Geometries

From the sf Cheatsheet

HW1 \(\rightarrow\) HW2

• Congrats on finishing HW1! You now know how to create geometries with sf and terra
• So now, what can you do with them?
• For example, we’d like to be able to say things like:
  • “The new lamppost cannot be placed at \((x, y)\), since there is already a building there!”
  • “There are \(N_1\) lampposts in County 1, and \(N_2\) lampposts in County 2”
  • “The average resident in Neighborhood A lives 2 km away from their nearest bus stop

First Things First: Loading and Saving

• Note how there were no data files in HW1 😱
• From HW2 onwards (and in your GIS life), we’ll:
  • Download from e.g. city Open Data Portals: geo data files, but also loading on-the-fly (this week)
  • Summarize/aggregate (this week and next week)
  • Visualize findings (“Mapping Libraries” unit)

Vector Formats

Shapefiles (.shp et al.)

A shape“file” is actually (at least) three separate files bundled together:

• Mandatory .shp: Containing feature geometries
• Mandatory .shx: Positional indices
• Mandatory .dbf: Data attributes
• Optional .prj: Coordinate reference system
• Optional .xml: Metadata

Shapefiles

Let’s see what’s inside the shapefile we first saw in Week 1, containing data on DC’s Census Tracts: Census Tracts in 2020

DC Census Tracts (with the Georgetown campus tract highlighted!) from OpenData.DC.gov

Shapefile Anatomy

From Rodrigue (2016)

GeoJSON / TopoJSON (.geojson)

• JavaScript Object Notation: General cross-platform format
• Useful when data is too complex for e.g. .csv
• TopoJSON = Memory-efficient GeoJSON
• Bonus: Inline preview on GitHub!

my_data.geojson

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [30, 20], [45, 40],
            [10, 40], [30, 20]
          ]
        ]
      },
      "properties": {
        "color": "green",
        "area": 3565747
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [15, 5], [40, 10],
            [10, 20], [5, 10], 
            [15, 5]
          ]
        ]
      },
      "properties": {
        "color": "red",
        "area": 3272386
      }
    }
  ]
}

GeoPackage (.gpkg)

• Open-source (non-proprietary) data format standard

Raster Formats

GeoTIFF (.tif)

• Based on TIFF format developed at NASA.

NetCDF (.nc4)

• Used in earth sciences, as format for data sources measured and distributed multiple times per day over large full-country or full-continent areas.

Coordinate Reference Systems (CRS)

• EPSG (European Petroleum Survey Group) Registry: Most common way to specify a CRS

  • For example, 4326 is the EPSG code for the WGS84 coordinate system

• PROJ: Rather than opaque numeric code like EPSG, uses plaintext “proj-strings” containing parameter info: datum, ellipsoid, projection, and units (e.g. meters). Example: PROJ4 code EPSG:4326 is represented as

  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

• WKT: Lengthy but human-readable descriptions

Geospatial Operations

Here we’ll look at Africa:

• The least “stretched” by different projections,
• Interesting geometric properties like countries-inside-of-countries (South Africa, Gambia), and
• Will be helpful for comparing/contrasting different ways of defining “neighbors”/neighborhoods

Getting the Geometries

Using rnaturalearth with mapview

France

library(rnaturalearth)
library(mapview)
france <- ne_countries(country = "France", scale = 50)
(france_map <- mapview(france, label = "geounit", legend = FALSE))

Centroid of France

library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

france_cent <- sf::st_centroid(france)

Warning: st_centroid assumes attributes are constant over geometries

france_map + mapview(france_cent, label = "Centroid", legend = FALSE)

One We Already Saw: Union

Computing the union of all geometries in the sf via sf::st_union()

library(leaflet.extras2)

Loading required package: leaflet

africa <- ne_countries(continent = "Africa", scale = 50)
africa_union <- sf::st_union(africa)
africa_map <- mapview(africa, label="geounit", legend=FALSE)
africa_union_map <- mapview(africa_union, label="st_union(africa)", legend=FALSE)
africa_map | africa_union_map

Helpful for Rasterizing: BBox

africa_bbox <- sf::st_bbox(africa)
africa_bbox_map <- mapview(africa_bbox, label="st_bbox(africa)", legend=FALSE)
africa_map | africa_bbox_map

Convex Hulls by Country

africa_countries_cvx <- sf::st_convex_hull(africa)
africa_countries_cvx_map <- mapview(africa_countries_cvx, label="geounit", legend=FALSE)
africa_map | africa_countries_cvx_map

Convex Hull of Continent

Use st_union() first:

africa_cvx <- africa |> st_union() |> st_convex_hull()
africa_cvx_map <- mapview(africa_cvx, label="geounit", legend=FALSE)
africa_map | africa_cvx_map

One We Already Saw: Centroids

Computing the centroid of all geometries in the sf via sf::st_centroid()

africa_cents <- sf::st_centroid(africa)

Warning: st_centroid assumes attributes are constant over geometries

africa_cents_map <- mapview(africa_cents, label="geounit", legend=FALSE)
africa_map | africa_cents_map

Pairs of Geometries

sa <- ne_countries(country = "South Africa", scale = 50)
lesotho <- ne_countries(country = "Lesotho", scale=50)
sf::st_relate(lesotho, sa)

although coordinates are longitude/latitude, st_relate assumes that they are
planar

     [,1]       
[1,] "FF2F1F212"

sf::st_covers(lesotho, sa)

Sparse geometry binary predicate list of length 1, where the predicate
was `covers'
 1: (empty)

sf::st_covered_by(lesotho, sa)

Sparse geometry binary predicate list of length 1, where the predicate
was `covered_by'
 1: (empty)

mapview(sa) + mapview(lesotho)

References

Rodrigue, Jean-Paul. 2016. The Geography of Transport Systems. Taylor & Francis.